Data inladen

KMIdata<-read.table(file = "~/GitHub/vespa_analyses/Input/KMIdata.txt", header=TRUE, sep="\t")

# Outlier vliegsnelheid v 10m/s weglaten
library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

KMIdata<-KMIdata%>%  filter(ID!=195)

library(ggplot2)
library(ggpubr)

Volgende datasets zijn gecreëerd:

KMIdata: dataset met alle vliegtijden en parameters (KMI)

Voor het model: Flight time ~ Distance gebruiken we de dataset vliegtijden per individu Omdat we hiermee de theoretische regel 1min=100m kunnen verifiëren. De imkers nemen hiervoor altijd kortste meting

Voor modellen met weerparameters, gewicht, urbanisatie: Hiervoor nemen we telkens de hele dataset, omdat elke meting van deze factoren afhangt.

Outlier van nest 29, Bait 1, Individu A weggelaten in deze dataset. (600m op 2min lijkt wel heel snel) Outlier van Melle ook weggelaten (meting op 2 km niet betrouwbaar) Meting in Excel script wel nog terug te vinden.

1.1) ForagingSpeed vs Temperature

Graphs

Van temperatuur ook categorische variabele maken (aanraden van Prof.Vangestel)

KMIdata$Temperature_cat<-cut(KMIdata$Temperature, c(8,13,17,21,25,29,33))

KMIdata$Temperature_cat

##   [1] (21,25] (21,25] (21,25] (21,25] (21,25] (21,25] (21,25] (21,25] (21,25]
##  [10] (21,25] (21,25] (21,25] (21,25] (21,25] (21,25] (21,25] (21,25] (17,21]
##  [19] (21,25] (17,21] (17,21] (17,21] (21,25] (21,25] (25,29] (25,29] <NA>   
##  [28] <NA>    <NA>    <NA>    <NA>    <NA>    (17,21] (25,29] (21,25] (17,21]
##  [37] (17,21] (25,29] (25,29] (25,29] (25,29] (25,29] (25,29] (25,29] (29,33]
##  [46] (25,29] (25,29] (25,29] (21,25] (21,25] (8,13]  (8,13]  (8,13]  (8,13] 
##  [55] (17,21] (8,13]  (8,13]  (8,13]  (8,13]  (8,13]  (8,13]  (17,21] (17,21]
##  [64] (17,21] (17,21] <NA>    (17,21] (17,21] (17,21] (17,21] (17,21] (17,21]
##  [73] (17,21] (17,21] (17,21] (17,21] (17,21] (17,21] (17,21] (17,21] (17,21]
##  [82] (25,29] (25,29] (25,29] (21,25] (21,25] (13,17] (13,17] (21,25] (21,25]
##  [91] (21,25] (21,25] (21,25] <NA>    <NA>    (17,21] (21,25] (21,25] <NA>   
## [100] <NA>    <NA>    <NA>    <NA>    <NA>    (13,17] (13,17] (13,17] (8,13] 
## [109] (8,13]  (8,13]  (25,29] (25,29] (13,17] (13,17] (17,21] (13,17] (13,17]
## [118] (13,17] (13,17] (13,17] (17,21] (17,21] (13,17] (13,17] (8,13]  (8,13] 
## [127] (8,13]  (13,17] (13,17] (17,21] (17,21] (13,17] (17,21] (17,21] <NA>   
## [136] (17,21] (13,17] (13,17] (13,17] (8,13]  (8,13]  (13,17] (13,17] (13,17]
## [145] (17,21] (17,21] (13,17] (17,21] (13,17] (13,17] (8,13]  (8,13]  (13,17]
## [154] (8,13]  (8,13]  (8,13]  (8,13]  (21,25] (21,25] (21,25] (17,21] (13,17]
## [163] (13,17] (13,17] (13,17] (13,17] (13,17] (13,17] (13,17] (13,17] (13,17]
## [172] (17,21] (17,21] (17,21] (17,21] (13,17] (17,21] (13,17] (13,17] (17,21]
## [181] (13,17] (8,13]  (13,17] (13,17] (13,17] (17,21] (17,21] (13,17] (13,17]
## [190] (8,13]  (8,13]  (8,13]  (8,13]  (8,13]  (8,13]  (8,13]  (8,13]  (8,13] 
## [199] (8,13]  (8,13]  (8,13]  (8,13]  (8,13]  (8,13]  (8,13]  (8,13]  (8,13] 
## [208] (17,21] (17,21] (17,21] (13,17] (13,17] (17,21] (17,21] (13,17] (13,17]
## [217] (17,21] (8,13]  (8,13]  (8,13]  (13,17] (13,17] (13,17] (13,17] (13,17]
## [226] (13,17] (8,13]  (8,13]  (8,13]  (8,13]  (8,13]  (8,13]  (8,13]  (8,13] 
## [235] (8,13]  (8,13]  (13,17] (13,17] (13,17] (13,17] (13,17] (8,13]  (8,13] 
## [244] (8,13]  (8,13]  (8,13]  (8,13]  (8,13]  (8,13]  (8,13]  (8,13]  (8,13] 
## [253] (13,17] (13,17] (17,21] (8,13]  (13,17] (13,17] (17,21] (17,21]
## Levels: (8,13] (13,17] (17,21] (21,25] (25,29] (29,33]

ggplot(data=subset(KMIdata, !is.na(Temperature_cat)), aes(x=Temperature_cat)) + geom_bar(fill="lightcoral")

Model Output

Temperature = continuous variable

Significant
Normality niet ok!

Temperature = discrete variable

Een aantal klassen significant, interpretatie?
Normality niet ok!

Poisson verdeling ook geprobeerd met glmer maar foutmelding (zie email)

library(lmerTest)

## Loading required package: lme4

## Loading required package: Matrix

## 
## Attaching package: 'lmerTest'

## The following object is masked from 'package:lme4':
## 
##     lmer

## The following object is masked from 'package:stats':
## 
##     step

library(lme4)

model_temp_cont<-lmer(Flighttime_min ~ Temperature_KMI + (1|Individualcode), offset=Distance, na.action=na.omit, data=KMIdata)
summary(model_temp_cont)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Flighttime_min ~ Temperature_KMI + (1 | Individualcode)
##    Data: KMIdata
##  Offset: Distance
## 
## REML criterion at convergence: 1898.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.8374 -0.2430 -0.0059  0.1753  3.8318 
## 
## Random effects:
##  Groups         Name        Variance  Std.Dev.
##  Individualcode (Intercept) 85390.122 292.216 
##  Residual                       3.288   1.813 
## Number of obs: 230, groups:  Individualcode, 91
## 
## Fixed effects:
##                  Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)     -416.0575    30.7203   91.0285 -13.543  < 2e-16 ***
## Temperature_KMI    0.3677     0.1250  138.0950   2.942  0.00382 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## Temprtr_KMI -0.075

res<-residuals(model_temp_cont)
shapiro.test(res)

## 
##  Shapiro-Wilk normality test
## 
## data:  res
## W = 0.83032, p-value = 3.923e-15

qqnorm(res)
qqline(res)

model_temp_cat<-lmer(Flighttime_min ~ Temperature_cat + (1|Individualcode), offset=Distance, data=KMIdata)
summary(model_temp_cat)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Flighttime_min ~ Temperature_cat + (1 | Individualcode)
##    Data: KMIdata
##  Offset: Distance
## 
## REML criterion at convergence: 1864.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.9276 -0.2865 -0.0056  0.1522  3.8469 
## 
## Random effects:
##  Groups         Name        Variance  Std.Dev.
##  Individualcode (Intercept) 86380.491 293.906 
##  Residual                       3.094   1.759 
## Number of obs: 230, groups:  Individualcode, 91
## 
## Fixed effects:
##                         Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)            -412.7019    30.9886   89.0891 -13.318  < 2e-16 ***
## Temperature_cat(13,17]    2.8778     0.6471  135.0087   4.447  1.8e-05 ***
## Temperature_cat(17,21]    2.6153     0.9093  135.0332   2.876  0.00468 ** 
## Temperature_cat(21,25]    2.8244     1.5162  135.0699   1.863  0.06466 .  
## Temperature_cat(25,29]    2.4914     1.6347  135.0732   1.524  0.12983    
## Temperature_cat(29,33]  102.5019   295.5400   89.0012   0.347  0.72954    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) T_(13, T_(17, T_(21, T_(25,
## Tmp_(13,17] -0.014                            
## Tmp_(17,21] -0.020  0.469                     
## Tmp_(21,25] -0.019  0.282  0.600              
## Tmp_(25,29] -0.018  0.261  0.556  0.734       
## Tmp_(29,33] -0.105  0.001  0.002  0.002  0.002

anova(model_temp_cat, ddf="Satterthwaite", type=3)

## Type III Analysis of Variance Table with Satterthwaite's method
##                 Sum Sq Mean Sq NumDF DenDF F value   Pr(>F)   
## Temperature_cat 64.298   12.86     5 122.3  4.1557 0.001603 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

res<-residuals(model_temp_cat)
shapiro.test(res)

## 
##  Shapiro-Wilk normality test
## 
## data:  res
## W = 0.83917, p-value = 1.037e-14

qqnorm(res)
qqline(res)

1.2) Flight time vs Distance + Temperature

Graph

Klopt dit? Kan dit ook met mixed model?

library(knitr)
library(kableExtra)

## Warning in !is.null(rmarkdown::metadata$output) && rmarkdown::metadata$output
## %in% : 'length(x) = 2 > 1' in coercion to 'logical(1)'

## 
## Attaching package: 'kableExtra'

## The following object is masked from 'package:dplyr':
## 
##     group_rows

library("plot3D")

x <- KMIdata$Distance
y <- KMIdata$Temperature_KMI
z <- KMIdata$Flighttime_min

fit <- lm(z ~ x + y, na.action=na.exclude)
x.pred <- seq(min(x[!is.na(x)]), max(x[!is.na(x)]), length.out = 20)
y.pred <- seq(min(y[!is.na(y)]), max(y[!is.na(y)]), length.out = 20)
xy <- expand.grid( x = x.pred, y = y.pred)
z.pred <- matrix(predict(fit, newdata = xy), 
                 nrow = 20, ncol = 20)
fitpoints <- predict(fit)

scatter3D(x, y, z, pch = 19, cex = 0.6, colvar=FALSE, col="dodgerblue3", theta = 210, phi = 10, bty="u", col.panel ="grey93", expand =0.4, col.grid = "white", xlab = "Distance", ylab = "Temperature", zlab = "Flight time", surf = list(x = x.pred, y = y.pred, z = z.pred,  
facets = TRUE, col=ramp.col(col = c("dodgerblue4", "seagreen2"), n = 100, alpha=0.8), fit = fitpoints, border="black"),main = "Flight time vs Distance + Temperature")

Met plotly

Had met deze link problemen met expand.grid https://stackoverflow.com/questions/38331198/add-regression-plane-to-3d-scatter-plot-in-plotly Dan maar grid dat van hierboven gebruikt, maar nu klopt er precies iets niet? Alle datapunten liggen boven het oppervlak.

library(plotly)

## 
## Attaching package: 'plotly'

## The following object is masked from 'package:ggplot2':
## 
##     last_plot

## The following object is masked from 'package:stats':
## 
##     filter

## The following object is masked from 'package:graphics':
## 
##     layout

library(tidyverse)

## ── Attaching packages
## ───────────────────────────────────────
## tidyverse 1.3.2 ──

## ✔ tibble  3.1.8     ✔ purrr   1.0.1
## ✔ tidyr   1.2.1     ✔ stringr 1.5.0
## ✔ readr   2.1.4     ✔ forcats 1.0.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ tidyr::expand()          masks Matrix::expand()
## ✖ plotly::filter()         masks dplyr::filter(), stats::filter()
## ✖ kableExtra::group_rows() masks dplyr::group_rows()
## ✖ dplyr::lag()             masks stats::lag()
## ✖ tidyr::pack()            masks Matrix::pack()
## ✖ tidyr::unpack()          masks Matrix::unpack()

KMIdataNoNA<-KMIdata %>%  filter(Flighttime_min!="NA")
KMIdataNoNA<-KMIdataNoNA %>%  filter(Temperature_KMI!="NA")

fig <- plot_ly(KMIdataNoNA, x = ~Distance, y = ~Temperature_KMI, z = ~Flighttime_min, size=1)
fig <- fig %>% add_markers()
fig <- fig %>% layout(scene = list(xaxis = list(title = 'Distance'),
                     yaxis = list(title = 'Temperature_KMI'),
                     zaxis = list(title = 'Flight time (min)')))
p2 <- add_trace(p = fig,
                z = z.pred,
                x = seq(2100, 0, by = -100),
                y = seq(0, 30, by = 10),
                type = "surface")
p2

Model Ouput

Normality niet ok!
Niet significant

model_tempdist_all<-lmer(Flighttime_min ~ Distance + Temperature_KMI + (1|NestID) + (1|NestID:BaitID) + (1|Individualcode), na.action=na.exclude, data=KMIdata)

## boundary (singular) fit: see help('isSingular')

summary(model_tempdist_all)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Flighttime_min ~ Distance + Temperature_KMI + (1 | NestID) +  
##     (1 | NestID:BaitID) + (1 | Individualcode)
##    Data: KMIdata
## 
## REML criterion at convergence: 1005.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.2425 -0.4262 -0.1822  0.1635  4.4765 
## 
## Random effects:
##  Groups         Name        Variance Std.Dev.
##  Individualcode (Intercept) 0.1443   0.3799  
##  NestID:BaitID  (Intercept) 1.8910   1.3751  
##  NestID         (Intercept) 0.0000   0.0000  
##  Residual                   3.2814   1.8115  
## Number of obs: 230, groups:  Individualcode, 91; NestID:BaitID, 68; NestID, 34
## 
## Fixed effects:
##                   Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)      3.5399788  0.9074782 66.9650223   3.901 0.000225 ***
## Distance         0.0056150  0.0008062 62.3323406   6.965 2.42e-09 ***
## Temperature_KMI -0.0686672  0.0445457 91.2883267  -1.541 0.126655    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Distnc
## Distance    -0.336       
## Temprtr_KMI -0.903 -0.010
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')

anova(model_tempdist_all, ddf="Satterthwaite", type=3)

## Type III Analysis of Variance Table with Satterthwaite's method
##                  Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
## Distance        159.173 159.173     1 62.332 48.5081 2.415e-09 ***
## Temperature_KMI   7.797   7.797     1 91.288  2.3762    0.1267    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

res<-residuals(model_tempdist_all)
shapiro.test(res)

## 
##  Shapiro-Wilk normality test
## 
## data:  res
## W = 0.7916, p-value < 2.2e-16

qqnorm(res)
qqline(res)

2) ForagingSpeed vs Cloudcoverage

Graphs

ggplot(KMIdata, aes(x=Cloudcoverage_KMI, y=ForagingSpeed)) + geom_point(color="lightblue") + geom_smooth(method="lm", formula =y ~ x, se=TRUE, fullrange=FALSE, level=0.95, color="darkblue") + ggtitle("All data KMI")+ theme(plot.title = element_text(hjust = 0.5, size=10))

Model Ouput

Licht significant
Normality niet ok!

model_cloud_all<-lmer(Flighttime_min ~ Cloudcoverage_KMI + (1|Individualcode), offset=Distance, na.action=na.omit, data=KMIdata)
summary(model_cloud_all)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Flighttime_min ~ Cloudcoverage_KMI + (1 | Individualcode)
##    Data: KMIdata
##  Offset: Distance
## 
## REML criterion at convergence: 1574.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.8876 -0.2653 -0.0075  0.1399  3.4796 
## 
## Random effects:
##  Groups         Name        Variance  Std.Dev.
##  Individualcode (Intercept) 79726.616 282.359 
##  Residual                       4.034   2.009 
## Number of obs: 191, groups:  Individualcode, 74
## 
## Fixed effects:
##                   Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)       -393.232     32.826   73.019 -11.979   <2e-16 ***
## Cloudcoverage_KMI    1.296      0.965  116.030   1.343    0.182    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## Cldcvrg_KMI -0.011

anova(model_cloud_all, ddf="Satterthwaite", type=3)

## Type III Analysis of Variance Table with Satterthwaite's method
##                   Sum Sq Mean Sq NumDF  DenDF F value Pr(>F)
## Cloudcoverage_KMI 7.2713  7.2713     1 116.03  1.8024  0.182

res<-residuals(model_cloud_all)
shapiro.test(res)

## 
##  Shapiro-Wilk normality test
## 
## data:  res
## W = 0.82333, p-value = 5.935e-14

qqnorm(res)
qqline(res)

3.1) ForagingSpeed vs Windspeed

Graphs

Telkens voor de modellen:

ForagingSpeed ~ Windspeed

ForagingSpeed ~ Windspeed²

plot19<-ggplot(KMIdata, aes(x=WindSpeed_KMI, y=ForagingSpeed)) + geom_point(color="khaki3") + geom_smooth(method="lm", formula =y ~ x, se=TRUE, fullrange=FALSE, level=0.95, color="olivedrab") + ggtitle("ForagingSpeed ~ Windspeed  KMI")+ theme(plot.title = element_text(hjust = 0.5, size=10))

plot22<-ggplot(KMIdata, aes(x=WindSpeed_KMI, y=ForagingSpeed)) + geom_point(color="khaki3") + geom_smooth(method="lm", formula =y ~ I(x^2), se=TRUE, fullrange=FALSE, level=0.95, color="olivedrab") + ggtitle("ForagingSpeed ~ Windspeed²  KMI")+ theme(plot.title = element_text(hjust = 0.5, size=10))


ggarrange(plot19, plot22 + rremove("x.text"), 
          labels = c("A", "B"),
          ncol = 2, nrow = 1)

Model Output

Beide modellen significant
Normality voor beiden niet ok!

model_wind_all<-lmer(Flighttime_min ~ WindSpeed_KMI + (1|Individualcode), offset=Distance, na.action=na.omit, data=KMIdata)
summary(model_wind_all)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Flighttime_min ~ WindSpeed_KMI + (1 | Individualcode)
##    Data: KMIdata
##  Offset: Distance
## 
## REML criterion at convergence: 1891.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.8936 -0.2627 -0.0058  0.1435  3.4706 
## 
## Random effects:
##  Groups         Name        Variance  Std.Dev.
##  Individualcode (Intercept) 85306.430 292.073 
##  Residual                       3.144   1.773 
## Number of obs: 230, groups:  Individualcode, 91
## 
## Fixed effects:
##                Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)   -411.5415    30.6235   90.0651 -13.439  < 2e-16 ***
## WindSpeed_KMI    0.6332     0.1610  138.0169   3.932 0.000133 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## WindSpd_KMI -0.019

anova(model_wind_all, ddf="Satterthwaite", type=3)

## Type III Analysis of Variance Table with Satterthwaite's method
##               Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
## WindSpeed_KMI 48.614  48.614     1 138.02  15.464 0.0001325 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

model_wind2_all<-lmer(Flighttime_min ~ I(WindSpeed_KMI^2) + (1|Individualcode), offset=Distance, na.action=na.omit, data=KMIdata)

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## unable to evaluate scaled gradient

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge: degenerate Hessian with 1 negative eigenvalues

summary(model_wind2_all)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Flighttime_min ~ I(WindSpeed_KMI^2) + (1 | Individualcode)
##    Data: KMIdata
##  Offset: Distance
## 
## REML criterion at convergence: 1902.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.1142 -0.2555 -0.0059  0.1360  3.6182 
## 
## Random effects:
##  Groups         Name        Variance  Std.Dev.
##  Individualcode (Intercept) 85336.122 292.123 
##  Residual                       3.302   1.817 
## Number of obs: 230, groups:  Individualcode, 91
## 
## Fixed effects:
##                      Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)        -410.12865   30.62483   90.01791 -13.392  < 2e-16 ***
## I(WindSpeed_KMI^2)    0.05677    0.01997  138.01433   2.844  0.00514 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## I(WS_KMI^2) -0.010
## optimizer (nloptwrap) convergence code: 0 (OK)
## unable to evaluate scaled gradient
## Model failed to converge: degenerate  Hessian with 1 negative eigenvalues

anova(model_wind2_all, ddf="Satterthwaite", type=3)

## Type III Analysis of Variance Table with Satterthwaite's method
##                    Sum Sq Mean Sq NumDF  DenDF F value  Pr(>F)   
## I(WindSpeed_KMI^2) 26.697  26.697     1 138.01  8.0858 0.00514 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

res<-residuals(model_wind_all)
shapiro.test(res)

## 
##  Shapiro-Wilk normality test
## 
## data:  res
## W = 0.81594, p-value = 8.702e-16

qqnorm(res)
qqline(res)

res<-residuals(model_wind2_all)
shapiro.test(res)

## 
##  Shapiro-Wilk normality test
## 
## data:  res
## W = 0.80097, p-value < 2.2e-16

qqnorm(res)
qqline(res)

3.2) ForagingSpeed vs Windspeed | afhankelijk van de windrichting

Hiervoor werd telkens voor elke meting nagegaan of de hoornaar met meewind (tailwind), tegenwind (upwind) of loodrechte wind (perpendicular) te maken had. Dit volgens de formules:

|𝜃flight −𝜃wind| ≤ 45 is tailwind

45 < |𝜃flight −𝜃wind|<135 is (quasi) perpendicular

|𝜃 flight −𝜃wind| ≥135 upwind

Graphs

TO DO: Lege vakjes Wind_flight NA maken

ggplot(data=subset(KMIdata, !is.na(Wind_flight)), aes(x= Wind_flight, col=Wind_flight, y=ForagingSpeed)) + geom_boxplot()

## Warning: Removed 26 rows containing non-finite values (`stat_boxplot()`).

ggplot(data=subset(KMIdata, !is.na(Wind_flight)), aes(x= WindSpeed_KMI, col=Wind_flight, y=ForagingSpeed)) + geom_point() + facet_grid(~Wind_flight) + geom_smooth(method="lm", formula = y ~ I(x^2), se=TRUE, fullrange=FALSE, level=0.95)

## Warning: Removed 30 rows containing non-finite values (`stat_smooth()`).

## Warning: Removed 30 rows containing missing values (`geom_point()`).

Model Output

Anova van Foragingspeed en de 3 windinvloeden

interpretatie?

windanova<-lmer(Flighttime_min ~ Wind_flight + (1|Individualcode), offset=Distance, na.action=na.exclude, data=KMIdata)

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## unable to evaluate scaled gradient

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge: degenerate Hessian with 1 negative eigenvalues

summary(windanova)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Flighttime_min ~ Wind_flight + (1 | Individualcode)
##    Data: KMIdata
##  Offset: Distance
## 
## REML criterion at convergence: 1932.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.8297 -0.2558 -0.0060  0.1122  4.1052 
## 
## Random effects:
##  Groups         Name        Variance  Std.Dev.
##  Individualcode (Intercept) 84703.912 291.039 
##  Residual                       3.441   1.855 
## Number of obs: 234, groups:  Individualcode, 94
## 
## Fixed effects:
##                          Estimate Std. Error      df t value Pr(>|t|)  
## (Intercept)               -337.11     168.03   91.99  -2.006   0.0478 *
## Wind_flightperpendicular   -72.48     170.78   91.99  -0.424   0.6723  
## Wind_flighttailwind        -71.34     170.78   92.00  -0.418   0.6771  
## Wind_flightupwind          -72.06     170.78   92.00  -0.422   0.6741  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Wnd_flghtpr Wnd_flghtt
## Wnd_flghtpr -0.984                       
## Wnd_flghttl -0.984  1.000                
## Wnd_flghtpw -0.984  1.000       1.000    
## optimizer (nloptwrap) convergence code: 0 (OK)
## unable to evaluate scaled gradient
## Model failed to converge: degenerate  Hessian with 1 negative eigenvalues

ForagingSpeed vs WindSpeed voor Tailwind dataset

Significant
Normality niet ok!

data_tailwind<-subset(KMIdata, Wind_flight == "tailwind")

model_tailwind<-lmer(Flighttime_min ~ WindSpeed_KMI + (1|Individualcode), offset=Distance, na.action=na.omit, data=data_tailwind)

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## unable to evaluate scaled gradient

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge: degenerate Hessian with 1 negative eigenvalues

summary(model_tailwind)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Flighttime_min ~ WindSpeed_KMI + (1 | Individualcode)
##    Data: data_tailwind
##  Offset: Distance
## 
## REML criterion at convergence: 372.3
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.13557 -0.15923 -0.00299  0.01975  2.13987 
## 
## Random effects:
##  Groups         Name        Variance  Std.Dev.
##  Individualcode (Intercept) 89183.261 298.64  
##  Residual                       1.323   1.15  
## Number of obs: 42, groups:  Individualcode, 22
## 
## Fixed effects:
##                Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)   -354.2999    63.6953   21.0337  -5.562  1.6e-05 ***
## WindSpeed_KMI    1.9695     0.5308   19.0056   3.710  0.00148 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## WindSpd_KMI -0.028
## optimizer (nloptwrap) convergence code: 0 (OK)
## unable to evaluate scaled gradient
## Model failed to converge: degenerate  Hessian with 1 negative eigenvalues

anova(model_tailwind, ddf="Satterthwaite", type=3)

## Type III Analysis of Variance Table with Satterthwaite's method
##               Sum Sq Mean Sq NumDF  DenDF F value   Pr(>F)   
## WindSpeed_KMI 18.206  18.206     1 19.006  13.765 0.001484 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

res<-residuals(model_tailwind)
shapiro.test(res)

## 
##  Shapiro-Wilk normality test
## 
## data:  res
## W = 0.80132, p-value = 4.725e-06

qqnorm(res)
qqline(res)

ForagingSpeed vs WindSpeed voor Perpendicular dataset

Significant
Normality niet ok!

data_perpendicular<-subset(KMIdata, Wind_flight == "perpendicular")


model_perpendicular<-lmer(Flighttime_min ~ WindSpeed_KMI + (1|Individualcode), offset=Distance, na.action=na.omit, data=data_perpendicular)
summary(model_perpendicular)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Flighttime_min ~ WindSpeed_KMI + (1 | Individualcode)
##    Data: data_perpendicular
##  Offset: Distance
## 
## REML criterion at convergence: 1086.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.6773 -0.3007 -0.0055  0.1599  3.8205 
## 
## Random effects:
##  Groups         Name        Variance  Std.Dev.
##  Individualcode (Intercept) 82549.958 287.315 
##  Residual                       2.418   1.555 
## Number of obs: 132, groups:  Individualcode, 55
## 
## Fixed effects:
##                Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)   -392.4157    38.7470   54.0285 -10.128 4.33e-14 ***
## WindSpeed_KMI    0.7233     0.1691   76.0084   4.277 5.46e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## WindSpd_KMI -0.016

anova(model_perpendicular, ddf="Satterthwaite", type=3)

## Type III Analysis of Variance Table with Satterthwaite's method
##               Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
## WindSpeed_KMI  44.22   44.22     1 76.008  18.289 5.456e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

res<-residuals(model_perpendicular)
shapiro.test(res)

## 
##  Shapiro-Wilk normality test
## 
## data:  res
## W = 0.84239, p-value = 1.443e-10

qqnorm(res)
qqline(res)

ForagingSpeed vs WindSpeed voor Upwind dataset

Niet significant
Normality niet ok!

data_upwind<-subset(KMIdata, Wind_flight == "upwind")


model_upwind<-lmer(Flighttime_min ~ WindSpeed_KMI + (1|Individualcode), offset=Distance, na.action=na.omit, data=data_upwind)
summary(model_upwind)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Flighttime_min ~ WindSpeed_KMI + (1 | Individualcode)
##    Data: data_upwind
##  Offset: Distance
## 
## REML criterion at convergence: 511.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.9621 -0.2480 -0.0058  0.1086  2.9613 
## 
## Random effects:
##  Groups         Name        Variance  Std.Dev.
##  Individualcode (Intercept) 89543.074 299.237 
##  Residual                       4.592   2.143 
## Number of obs: 56, groups:  Individualcode, 27
## 
## Fixed effects:
##                Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)   -412.9125    57.6072   26.0315  -7.168 1.29e-07 ***
## WindSpeed_KMI   -0.2379     0.4400   28.0061  -0.541    0.593    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## WindSpd_KMI -0.025

anova(model_upwind, ddf="Satterthwaite", type=3)

## Type III Analysis of Variance Table with Satterthwaite's method
##               Sum Sq Mean Sq NumDF  DenDF F value Pr(>F)
## WindSpeed_KMI 1.3425  1.3425     1 28.006  0.2924  0.593

res<-residuals(model_upwind)
shapiro.test(res)

## 
##  Shapiro-Wilk normality test
## 
## data:  res
## W = 0.80469, p-value = 3.644e-07

qqnorm(res)
qqline(res)

Model Selection

fullmodel<-lmer(Flighttime_min ~ Traject100m + WindSpeed_KMI + Temperature_KMI + Cloudcoverage_KMI + Weight_ind + Wind_flight + (1|Individualcode), offset=Distance, na.action=na.omit, data=KMIdata)
summary(fullmodel)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Flighttime_min ~ Traject100m + WindSpeed_KMI + Temperature_KMI +  
##     Cloudcoverage_KMI + Weight_ind + Wind_flight + (1 | Individualcode)
##    Data: KMIdata
##  Offset: Distance
## 
## REML criterion at convergence: 410.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.3917 -0.3001 -0.0556  0.2722  3.3894 
## 
## Random effects:
##  Groups         Name        Variance  Std.Dev.
##  Individualcode (Intercept) 53883.948 232.129 
##  Residual                       3.547   1.883 
## Number of obs: 71, groups:  Individualcode, 16
## 
## Fixed effects:
##                      Estimate Std. Error        df t value Pr(>|t|)   
## (Intercept)         -550.3334   308.0425   13.0598  -1.787  0.09723 . 
## Traject100m         1673.6657   550.0121   12.9981   3.043  0.00943 **
## WindSpeed_KMI         -0.1017     1.6066   50.0768  -0.063  0.94978   
## Temperature_KMI        0.5069     0.6809   50.0637   0.744  0.46009   
## Cloudcoverage_KMI     -1.2680     2.4523   50.0386  -0.517  0.60739   
## Weight_ind          -479.3634   861.3756   13.0048  -0.557  0.58731   
## Wind_flighttailwind    0.9777     1.3264   50.0094   0.737  0.46447   
## Wind_flightupwind      0.4210     1.0123   49.9981   0.416  0.67929   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Trj100 WS_KMI Tm_KMI Cl_KMI Wght_n Wnd_flghtt
## Traject100m -0.352                                              
## WindSpd_KMI  0.050 -0.019                                       
## Temprtr_KMI -0.052  0.017 -0.930                                
## Cldcvrg_KMI -0.031  0.016 -0.591  0.558                         
## Weight_ind  -0.874 -0.105 -0.025  0.025  0.012                  
## Wnd_flghttl -0.003  0.006 -0.220  0.001  0.303  0.001           
## Wnd_flghtpw -0.010  0.006 -0.188  0.158  0.296  0.003  0.238

model1<-lmer(Flighttime_min ~ Traject100m + Temperature_KMI + Cloudcoverage_KMI + Weight_ind + Wind_flight + (1|Individualcode), offset=Distance, na.action=na.omit, data=KMIdata)
summary(model1)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Flighttime_min ~ Traject100m + Temperature_KMI + Cloudcoverage_KMI +  
##     Weight_ind + Wind_flight + (1 | Individualcode)
##    Data: KMIdata
##  Offset: Distance
## 
## REML criterion at convergence: 413.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.4009 -0.3049 -0.0583  0.2752  3.4189 
## 
## Random effects:
##  Groups         Name        Variance  Std.Dev.
##  Individualcode (Intercept) 53842.509 232.040 
##  Residual                       3.479   1.865 
## Number of obs: 71, groups:  Individualcode, 16
## 
## Fixed effects:
##                      Estimate Std. Error        df t value Pr(>|t|)   
## (Intercept)         -549.3616   307.5411   13.0053  -1.786  0.09737 . 
## Traject100m         1673.0213   549.7055   12.9996   3.043  0.00942 **
## Temperature_KMI        0.4667     0.2481   51.0125   1.881  0.06569 . 
## Cloudcoverage_KMI     -1.3589     1.9593   51.0134  -0.694  0.49111   
## Weight_ind          -480.7074   860.7830   12.9992  -0.558  0.58603   
## Wind_flighttailwind    0.9596     1.2813   51.0107   0.749  0.45734   
## Wind_flightupwind      0.4090     0.9845   51.0036   0.415  0.67961   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Trj100 Tm_KMI Cl_KMI Wght_n Wnd_flghtt
## Traject100m -0.352                                       
## Temprtr_KMI -0.016 -0.001                                
## Cldcvrg_KMI -0.002  0.006  0.029                         
## Weight_ind  -0.874 -0.106  0.006 -0.004                  
## Wnd_flghttl  0.008  0.002 -0.568  0.220 -0.004           
## Wnd_flghtpw  0.000  0.003 -0.046  0.233 -0.002  0.206

model2<-lmer(Flighttime_min ~ Traject100m + Temperature_KMI + Cloudcoverage_KMI + Weight_ind + (1|Individualcode), offset=Distance, na.action=na.omit, data=KMIdata)
summary(model2)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Flighttime_min ~ Traject100m + Temperature_KMI + Cloudcoverage_KMI +  
##     Weight_ind + (1 | Individualcode)
##    Data: KMIdata
##  Offset: Distance
## 
## REML criterion at convergence: 418
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.5700 -0.2799 -0.0505  0.2098  3.4708 
## 
## Random effects:
##  Groups         Name        Variance Std.Dev.
##  Individualcode (Intercept) 53760.67 231.863 
##  Residual                       3.39   1.841 
## Number of obs: 71, groups:  Individualcode, 16
## 
## Fixed effects:
##                    Estimate Std. Error        df t value Pr(>|t|)   
## (Intercept)       -551.1020   307.2951   13.0051  -1.793  0.09618 . 
## Traject100m       1671.8808   549.2846   12.9996   3.044  0.00941 **
## Temperature_KMI      0.5674     0.2009   53.0059   2.825  0.00666 **
## Cloudcoverage_KMI   -1.7809     1.8501   53.0103  -0.963  0.34013   
## Weight_ind        -477.7494   860.1200   13.0005  -0.555  0.58802   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Trj100 Tm_KMI Cl_KMI
## Traject100m -0.352                     
## Temprtr_KMI -0.014  0.000              
## Cldcvrg_KMI -0.004  0.006  0.179       
## Weight_ind  -0.874 -0.106  0.004 -0.003

model3<-lmer(Flighttime_min ~ Traject100m + Temperature_KMI + Cloudcoverage_KMI + (1|Individualcode), offset=Distance, na.action=na.omit, data=KMIdata)
summary(model3)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Flighttime_min ~ Traject100m + Temperature_KMI + Cloudcoverage_KMI +  
##     (1 | Individualcode)
##    Data: KMIdata
##  Offset: Distance
## 
## REML criterion at convergence: 1536.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.5917 -0.2547 -0.0054  0.1572  3.2093 
## 
## Random effects:
##  Groups         Name        Variance  Std.Dev.
##  Individualcode (Intercept) 63054.182 251.106 
##  Residual                       3.757   1.938 
## Number of obs: 191, groups:  Individualcode, 74
## 
## Fixed effects:
##                    Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)       -559.5955    45.6541   72.7061 -12.257  < 2e-16 ***
## Traject100m       1018.4148   225.8033   72.0127   4.510 2.46e-05 ***
## Temperature_KMI      0.4846     0.1545  115.1085   3.138  0.00216 ** 
## Cloudcoverage_KMI    1.0641     0.9342  115.0364   1.139  0.25704    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Trj100 Tm_KMI
## Traject100m -0.766              
## Temprtr_KMI -0.070  0.012       
## Cldcvrg_KMI -0.002 -0.001 -0.079

model4<-lmer(Flighttime_min ~ Traject100m + Temperature_KMI  + (1|Individualcode), offset=Distance, na.action=na.omit, data=KMIdata)
summary(model4)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Flighttime_min ~ Traject100m + Temperature_KMI + (1 | Individualcode)
##    Data: KMIdata
##  Offset: Distance
## 
## REML criterion at convergence: 1867.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.8341 -0.2430 -0.0060  0.1773  3.8263 
## 
## Random effects:
##  Groups         Name        Variance  Std.Dev.
##  Individualcode (Intercept) 70771.693 266.029 
##  Residual                       3.288   1.813 
## Number of obs: 230, groups:  Individualcode, 91
## 
## Fixed effects:
##                  Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)     -559.5730    42.8324   89.6604 -13.064  < 2e-16 ***
## Traject100m     1009.9184   228.1899   89.0088   4.426 2.72e-05 ***
## Temperature_KMI    0.3728     0.1250  138.0877   2.983  0.00338 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Trj100
## Traject100m -0.757       
## Temprtr_KMI -0.061  0.009

library(cAIC4)

## Loading required package: stats4

## Loading required package: nlme

## 
## Attaching package: 'nlme'

## The following object is masked from 'package:lme4':
## 
##     lmList

## The following object is masked from 'package:dplyr':
## 
##     collapse

library(MASS)

## 
## Attaching package: 'MASS'

## The following object is masked from 'package:plotly':
## 
##     select

## The following object is masked from 'package:dplyr':
## 
##     select

library(AICcmodavg)

## 
## Attaching package: 'AICcmodavg'

## The following object is masked from 'package:lme4':
## 
##     checkConv

AIC(fullmodel)

## [1] 430.4662

AIC(model1)

## [1] 431.2467

AIC(model2)

## [1] 431.9586

AIC(model3)

## [1] 1548.611

AIC(model4)

## [1] 1877.757

Multicollinearity check

library(ggcorrplot)
library(car)

## Loading required package: carData

## 
## Attaching package: 'car'

## The following object is masked from 'package:purrr':
## 
##     some

## The following object is masked from 'package:dplyr':
## 
##     recode

datacor<-KMIdata[ , c("Temperature_KMI", "Cloudcoverage_KMI", "WindSpeed_KMI", "Weight_ind", "Traject100m", "NestID")]
source("~/GitHub/vespa_analyses/Input/HighstatLibV10.R") 
corvif(datacor)

## 
## 
## Variance inflation factors
## 
##                       GVIF
## Temperature_KMI   3.538307
## Cloudcoverage_KMI 1.333698
## WindSpeed_KMI     3.811471
## Weight_ind        2.981154
## Traject100m       1.435785
## NestID            5.380406

cormat <- round(cor(datacor, use = "pairwise.complete.obs"), 2)
ggcorrplot(cormat, lab= TRUE, type = "lower", ggtheme = ggplot2::theme_gray,
   colors = c("#6D9EC1", "white", "#E46726"))

Masterthesis script,weerparameters KMI

Heleen Van Ransbeeck

2023-02-27

Data inladen

1.1) ForagingSpeed vs Temperature

Graphs

Model Output

1.2) Flight time vs Distance + Temperature

Graph

Model Ouput

2) ForagingSpeed vs Cloudcoverage

Graphs

Model Ouput

3.1) ForagingSpeed vs Windspeed

Graphs

Model Output

3.2) ForagingSpeed vs Windspeed | afhankelijk van de windrichting

Graphs

Model Output

Anova van Foragingspeed en de 3 windinvloeden

ForagingSpeed vs WindSpeed voor Tailwind dataset

ForagingSpeed vs WindSpeed voor Perpendicular dataset

ForagingSpeed vs WindSpeed voor Upwind dataset

Model Selection

Multicollinearity check